********************************************************************************
// Robustness checks
********************************************************************************

clear all
set more off 
cap log close
log using "BR_EntryThreshold_PCSA_Robust.log", replace

capture program drop brentry_bioprobit_neg
qui do "brentry_bioprobit_negbin_v3.do"

putexcel set "BR_EntryThreshold_Robust_V5.xlsx", modify sheet("entry")
putexcel A1 = "UCC Licensing"
putexcel B1 = "NY and NC"
putexcel C1 = "Year 2014"
putexcel D1 = "Year 2016"
putexcel E1 = "Hosp w/ ED"
putexcel F1 = "Drop Texas"
putexcel G1 = "Urgent in name"
putexcel H1 = "HSA"
putexcel I1 = "All UCCs"

constraint drop
constrain 1 tot_pop = 1
constrain 2 tot_pop2 = 1

// Without UCC licensing

use "PCSALevelData_v3.dta", clear
drop if ucc_regulation==1
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo clear
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "A"

// Exclude NY and NC

use "PCSALevelData_v3.dta", clear
drop if inlist(state, "36", "37")
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "B"

estout *, drop(tot_pop*) cells(b se) label mlabels("UCC regulation" "NY and NC")

// Alternative years

forvalues y = 2014(2)2016{

	use "ZCTALevelRobust_`y'.dta", clear
	collapse (sum) tot_pop n_urgentcare n_hospitals n_hospaffucc_geo n_hospaffucc_geo2 any_emergency hisp black other white og_* (firstnm) state (max) ucc_regulation, by(pcsa)
	tempfile base 
	save `base'

	use "ZCTALevelProcessed_R2_input.dta", clear
	collapse (mean) rural income_pc cms_wage_index median_value median_gross_rent con_intensity pop_growth rpl_themes [w=tot_pop], by(pcsa)
	merge 1:1 pcsa using `base', nogen
	save "PCSALevelData_Raw_`y'.dta", replace 

	use "PCSALevelData_Raw_`y'.dta", clear
	merge 1:1 pcsa using "market_definition.dta", keep(3) nogen
	gen le_highschool = og_less_highschool + og_highschool
	gen gte_highschool = og_highschool + og_some_college + og_bachelor
	gen female = og_female 
	gen age_65 = og_age_65 
	gen uninsured = og_uninsured
	ren hisp hispanic 
	ren black nonhisp_black 
	ren other nonhisp_other
	ren white nonhisp_white
	replace tot_pop = tot_pop*10000*1000
	local vars "hispanic nonhisp_white nonhisp_black nonhisp_other female age_65 uninsured le_highschool gte_highschool"
	foreach v of local vars{
		 replace `v' = `v'/tot_pop
	}
	replace tot_pop = tot_pop/10000/1000
	drop if income_pc==.

	gen cat_ucc = n_urgentcare
	replace cat_ucc = 3 if cat_ucc>3
	gen cat_hosp = n_hospitals
	replace cat_hosp = 1 if cat_hosp>1
	gen cat_hosp2 = n_hospitals
	replace cat_hosp2 = 0 if cat_hosp2<=1
	replace cat_hosp2 = 1 if cat_hosp2>=2
	gen cat_aucc = n_hospaffucc_geo
	replace cat_aucc = 1 if cat_aucc>1
	gen n_ucc_aucc = n_urgentcare + n_hospaffucc_geo
	gen cat_both = n_ucc_aucc
	replace cat_both = 4 if cat_both>4
	gen tot_pop2 = tot_pop
	gen tot_pop3 = tot_pop
	gen any_hosp = (n_hospitals>0)
	gen any_aucc = (n_hospaffucc_geo>0)
	gen any_ucc = (n_urgentcare>0) 

	gen median_income_pc = .
	qui sum income_pc, detail
	replace median_income_pc = r(p50)
	gen high_income = (income_pc>=median_income_pc)
	drop median_income_pc

	gen median_svi = .
	qui sum rpl_themes, detail
	replace median_svi = r(p50)
	gen high_svi = (rpl_themes>=median_svi)
	drop median_svi

	gen median_uninsured = .
	qui sum uninsured, detail
	replace median_uninsured = r(p50)
	gen high_uninsured = (uninsured>=median_uninsured)
	drop median_uninsured

	ren n_hospitals og_n_hospitals
	gen n_hospitals = og_n_hospitals
	replace n_hospitals = cat_hosp2

	ren n_hospaffucc_geo og_n_hospaffucc_geo
	gen n_hospaffucc_geo = og_n_hospaffucc_geo
	replace n_hospaffucc_geo = cat_aucc

	ren n_urgentcare og_n_urgentcare 
	gen n_urgentcare = og_n_urgentcare
	replace n_urgentcare = cat_ucc

	label var rural "Rural"
	label var n_hospitals "Additional hospital presence"
	label var income_pc "Income per capita"
	label var hispanic "Hispanic"
	label var nonhisp_white "White"
	label var nonhisp_black "Black"
	label var gte_highschool "High school or more"
	label var age_65 "Age 65 or more"
	label var uninsured "Uninsured"
	label var n_urgentcare "Number of UCCs"
	label var n_hospaffucc_geo "Number of AUCCs"
	label var cms_wage_index "CMS wage index"
	label var con_intensity "CON laws" 

	save "PCSALevelData_v3_`y'.dta", replace 

}

use "PCSALevelData_v3_2014.dta", clear
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo clear
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "C"

use "PCSALevelData_v3_2016.dta", clear
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2) technique(bfgs)
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "D"

estout *, drop(tot_pop*) cells(b se) label mlabels("2014" "2016")

// Drop hospitals without emergency care

use "PCSALevelData_v3.dta", clear
replace cat_hosp2 = any_emergency
replace cat_hosp2 = 0 if cat_hosp2<=1
replace cat_hosp2 = 1 if cat_hosp2>=2
replace n_hospitals = cat_hosp2
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo clear
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "E"
estout *, drop(tot_pop*) cells("b(fmt(1)) se(fmt(1))") label

// Drop texas

use "PCSALevelData_v3.dta", clear
drop if state=="48"
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo clear
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "F"
estout *, drop(tot_pop*) cells("b(fmt(1)) se(fmt(1))") label

// HSA market definition

use "HSALevelDataProcessed_R2.dta", clear
drop cat_ucc cat_hosp
gen cat_ucc = n_urgentcare
replace cat_ucc = 3 if cat_ucc>3
gen cat_hosp = n_hospitals
replace cat_hosp = 1 if cat_hosp>1
gen cat_hosp2 = n_hospitals
replace cat_hosp2 = 0 if cat_hosp2<=1
replace cat_hosp2 = 1 if cat_hosp2>=2
gen cat_aucc = n_hospaffucc_geo
replace cat_aucc = 1 if cat_aucc>1
gen n_ucc_aucc = n_urgentcare + n_hospaffucc_geo
gen cat_both = n_ucc_aucc
replace cat_both = 4 if cat_both>4
gen any_aucc = (n_hospaffucc_geo>0)
gen any_ucc = (n_urgentcare>0) 
ren n_hospitals og_n_hospitals
gen n_hospitals = og_n_hospitals
replace n_hospitals = cat_hosp2
ren n_hospaffucc_geo og_n_hospaffucc_geo
gen n_hospaffucc_geo = og_n_hospaffucc_geo
replace n_hospaffucc_geo = cat_aucc
ren n_urgentcare og_n_urgentcare 
gen n_urgentcare = og_n_urgentcare
replace n_urgentcare = cat_ucc

ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "H"

// Including all UCCs 

use "PCSALevelData_v3.dta", clear
replace n_urgentcare = og_n_urgentcare + og_n_hospaffucc_geo 
replace cat_ucc = n_urgentcare 
replace cat_ucc = 3 if cat_ucc>3
replace n_urgentcare = cat_ucc
ml clear
ml model lf brentry_bioprobit_neg (hosp_s:cat_hosp2 = tot_pop, nocons) (hosp_v:cat_hosp2 = rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (hosp_f:cat_hosp2 = cms_wage_index con_intensity, nocons) /hosp_a1 /hosp_g1 (ucc_s:cat_ucc = tot_pop2, nocons) (ucc_v:cat_ucc = n_hospitals rural income_pc hispanic nonhisp_black gte_highschool age_65 uninsured, nocons) (ucc_f:cat_ucc = cms_wage_index, nocons) /ucc_a1 /ucc_a2 /ucc_a3 /ucc_g1 /ucc_g2 /ucc_g3 /r, constraints(1 2)
eststo clear
eststo: ml max, difficult iterate(200) tolerance(1e-4) ltolerance(1e-5)
qui do "BR_EntryThreshold_Bi.do" 10000 10000 "I"
estout *, drop(tot_pop*) cells("b(fmt(1)) se(fmt(1))") label

log close